AD-A275  648 

IIHIIIIHII 


© 


NRL  -  Contract  •  Reporting  information  for  the  months  of  August* 
September*  October,  November*  December  1993  and  January  1994. 


LOW  VOLTAGE  ELECTRON  BEAM  LITHOGRAPHY 
Aaron  Baum  and  R.  Fabian  Pease 


DTIC 


ELECTE 
FEB  8  1994 


In  order  to  develop  a  low-energy-spread,  high  brightness  electron 
source  using  negative  electron  affinity  technology,  it  is  necessary  to  survey 
the  effects  of  cathode  structure  and  activation  on  energy  spread,  lateral 
velocity  distribution,  peak  current,  and  lifetime.  Most  research  in  NEAPC 
technology  has  gone  into  detector  applications,  to  improve  photoyield  at  low 
light  intensity  and  high  wavelength.  Intevac  is  a  leader  in  this  area.  The 
scientific  understanding  and  technical  expertise  developed  in  this  effort  are 
important  to  achieving  our  goal  of  an  electron  source  optimized  for  low- 
energy  electron  applications;  however,  it  will  be  necessary  to  investigate  other 
cathode  properties  (peak  brightness,  energy  spread)  under  a  substantially 
different  operating  regime  (high  light  intensity,  small  emission  area,  high 
extraction  field). 


To  this  end  it  is  important  to  examine  the  operation  of  the 
photocathode  in  detail;  the  best  way  to  describe  this  is  by  tracking  the  life  of  an 
electron,  illustrated  in  Figure  1.  We  start  with  an  electron  in  the  valence 
band  (step  1),  which  is  excited  by  an  incoming  photon  into  the  conduction 
band  (step  2).  The  electron  then  relaxes,  by  optical  phonon  scattering  (which 
has  a  mean  free  path  of  about  300A),  into  h  e  conduction  band  minimum, 
where  the  electron  has  a  fairly  long  lifetime,  which  leads  to  a  long  diffusion 

length  (a  few  pm  in  a  good  photocathode).  Thus  electrons  excited  within  a 

few  pm  of  the  surface  have  a  high  probability  of  reaching  the  surface.  It 
should  be  noted  here  that  Intevac  has  developed  the  technology  of  a  glass- 
bonded  thin  film  photocathode  which  can  be  used  in  transmission  mode 
(light  entering  from  beneath  the  surface).  Besides  offering  ruggedness  and 
simplicity  of  electron  gun  design,  transmission  mode  NEAPCs  offer  lower 
energy  spread  than  the  more  commonly  used  reflection-mode  NEAPCs.  In 
reflection  mode,  there  are  many  more  "hot"  electrons  —  electrons  which 
haven’t  yet  thermalized  to  the  CMB  —  at  the  surface;  in  reflection  mode 
NEAPCs,  many  of  these  "hot"  electrons  escape,  broadening  the  energy  spread. 


When  the  electrons  reach  the  surface,  they  encounter  the  band-bending 
region,  where  they  are  accelerated  toward  the  surface.  Since  the  electrons  are 
"hot"  in  this  region,  they  may  interact  with  optical  phonons  here  and  lose  (or 
gain)  energy  (step  3).  Thus  their  energy  spread  increases.  Electrons  with 
energy  above  the  vacuum  level  can  then  escape.  The  work-function- 
lowering  activation  layer,  by  determining  the  vacuum  level,  acts  as  an  energy 
filter,  blocking  the  lower-energy  electrons  from  escaping. 
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There  is  another  phenomenon  at  the  surface  that  plays  a  potentially 
important  role  in  improving  brightness  of  NEAFC  sources.  Because  the 
electrons  have  appreciably  lower  effective  mass  in  the  semiconductor  than  in 
vacuum  (a  ratio  of  1:15  for  GaAs),  their  semicondcutor  quantum-mechanical 
wavelength  is  much  shorter.  When  the  electrons  are  emitted  into  vacuum, 

quantum-mechanical  refraction  takes  place,  and  the  electrons  are  focused  into 
a  narrow  forward  cone.  This  is  an  excellent  property  for  an  electron  source. 

The  activation  layer  also  plays  a  role  in  determining  emission 
characterisitics.  In  the  layer,  electrons  may  scatter  elastically,  which  would 
degrade  their  angular  distribution,  or  inelastically,  which  would  degrade  both 
their  angular  and  energy  distributions.  This  consideration  favors  a  thinner 
activation  layer. 

To  optimize  the  NEAPC  source  for  low-energy  eletron  beam  work,  we 
are  examining  the  effect  of  activation  layer  thickness  on  both  the  lateral  and 
total  energy  spreads.  One  structure  being  examined  is  a  GaAsP  cathode 
activated  with  only  a  monolayer  of  Cs.  This  type  of  NEAPC  has  not  been  fully 
investigated  as  it  has  lower  quantum  efficiency  and  less  sensitivity  to  long 
wavelengths  than  currently  used  cathodes.  However,  the  thin  activating 
layer  should  minimize  scattering,  improving  brightness,  and  an  extremely 
low  energy  spread  should  be  obtainable  by  adjusting  the  bandgap  and  using 
the  vacuum  level  as  an  energy  filter.  Other  possibilities  include  activation 
with  F  instead  of  O,  which  may  provide  thinner  activation  layers,  and 
increasing  the  doping  near  the  surface.  The  latter  modification  would 
decrease  the  width  of  the  depletion  region  and  thus  improve  the  energy 
spread  of  emitted  electrons  as  they  w;’!  be  less  likely  to  scatter  there. 

Figure  2  shows  a  possible  optimized  photocathode  structure.  Like  the 
typical  Intevac  cathodes,  this  structure  uses  AlGaAs  to  block  electrons  from 
diffusing  too  far  from  the  surface;  this  design  also  incorporates  a  graded 
bandgap  to  accelerate  the  electrons  toward  the  surface.  This  action  improves 
the  cathode's  efficiency  and  decreases  the  cathode's  response  time  (which  is  in 
the  tens  of  picoseconds  even  without  this  improvement).  Also,  since  a  small 
source  size  is  desireable  from  an  electron  optics  standpoint,  if  the  laser  is 
focused  to  a  diffraction-limited  spot  in  the  active  region,  the  graded  bandgap 
will  reduce  the  diffusive  spreading  of  the  electrons  before  emission.  To 
further  improve  the  source  size  and  response  time  of  the  cathode,  it  is 
designed  to  be  thinner  than  an  ordinary  cathode,  perhaps  less  than  1  micron. 
The  most  important  modifications,  however,  are  near  the  surface.  The 
depletion  region  has  been  reduced  by  "spike  doping,"  —  the  last  100  A  of  the 
cathode  are  doped  as  heavily  as  possible.  Furthermore,  a  thin  Cs-only 
activation  layer  is  employed.  These  modifications  should  greatly  reduce 
scattering  in  these  two  regions,  improving  energy  spread  and  brightness. 


Currently  we  are  examining  the  effect  of  activation  on  cathode 
brightness.  The  experiments  are  carried  out  using  photocathodes  made  using 
Intevac's  proprietary  techniques  sealed  in  tubes  at  UHV.  Our  first  series  of 
measurements  consists  of  measuring  the  lateral  energy  distributions  of 
electrons  emitted  from  various  photocathodes  by  monitoring  the  current 
intercepted  by  a  knife  edge  that  is  transported  across  the  beam.  The  tube 
design  is  shown  in  Figure  2.  This  apparatus  affords  the  highest  resolution 
measurements  of  this  type  ever  made  on  negative  electron  affinity  devices. 
Initial  measurements  have  been  made  on  photocathodes  activated  with  the 
standard  Intevac  nightvision  activation,  and  show  that,  as  expected,  the 
electrons  have  very  low  lateral  velocities,  corresponding  to  approximately 
40meV  on  average.  Tubes  with  GaAsP  cathodes  and  Cs-only  activations  are 
being  fabricated. 

Future  experiments  will  involve  measuring  total  energy  spread, 
lifetime,  and  peak  brightness  of  a  variety  of  cathode  structures. 
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FIRST  DRAFT _ 

Current  Status  of  Computer  Models  for  Charged  Particle 

Systems 


Brent  Boyer  end  R.  F.  W.  Pease 


Although  there  are  many  analytical  methods  available  for  analyzing  charged  particle  systems,  in 
general  they  apply  to  geometrically  simple  lens  systems.  As  a  result,  computer  modeling  has 
emerged  as  a  powerful  tool  for  designing  systems  which  will  extract,  focus,  and  analyze  charged 
particle  beams.  Personal  computers  are  now  sufficiently  powerful  and  economical  tha,  any 
researcher  can  have  access  to  the  necessary  computer  resources.  This  has  caused  computer 
modeling  to  become  very  popular.  Here  we  summarize  the  main  features  of  computer  modeling 
for  charged  particle  systems.  The  topics  covered  include  numerical  techniques,  source  modeling, 
lens  design,  deflector  modeling,  tolerances,  space  charge  effects,  and  beam-target  interactions. 


Introduction 


Although  the  earlier  textbooks  on  electron  (and  ion)  optics  antedated  the  widespread  use  of  computers,  die 
use  of  computers  for  solving  electron  optical  problems  goes  back  at  least  four  decades.  During  die  early  fifties 
Liebmann  and  his  colleagues  at  the  AEI  Research  Centre  in  the  UK  employed  a  combination  of  resistance  network 
analog  computation  and  digital  integration  to  solve  for  fields  and  published  a  series  of  papers  that  became  the  Trible' 
for  those  designing  axially  symmetric  electrostatic  and  magnetic  lenses  [I].  Munro  in  1970  described  a  series  of 
programs  for  doing  the  same  calculations  on  a  'modem'  digital  computer  far  more  rapidly  and  conveniently  [2).  His 
results  bore  out  those  of  the  AEI  team  almost  exactly. 


For  analyzing  such  lenses  the  first  task  is  solve  the  Laplace  equation  according  to  the  boundary  conditions 
set  by  the  magnetic  polepieces  and  their  magneto-static  potentials,  and  by  the  electrodes  and  their  electric  potentials. 
For  axial  symmetry  and  with  no  polepieces  or  electrodes  on  axis  (the  usual  arrangement)  the  properties  of  the 
solutions  to  the  Laplace  equation  are  such  that  it  is  possible  to  determine  both  first  order  properties  (focal  lengths, 
principal  plane  positions  and  chromatic  aberration  coefficient)  and  third  order  properties  (e.g.  spherical  aberration)  of 
the  resulting  lenses  from  only  the  on-axis  distributions  of  electric  and  magnetostatic  potentials.  V(0,z)  and  4(0,z). 

A  worked  example  for  a  magnetic  lens  is  shown  below  in  the  section  on  Lens  Modeling  using  a  series  of  programs 
written  by  Munro.  Munro's  programs  can  also  be  used  to  determine  the  magnetostatic  potential  when  the  material  in 
the  magnetic  circuit  approaches  saturation. 

The  limitation  of  the  above  treatment  is  that  it  does  not  take  into  account  higher  order  aberrations  (5th  order 
and  above)  which  may  be  significant  for  certain  cases.  To  do  a  precise  simulation  of  charged  particle  optics  requires 
determination  of  the  fields  everywhere,  not  just  on-axis,  and  the  particle  trajectories  are  determined  by  nnmmnii 
integration  of  the  corresponding  equation  of  motion.  The  techniques  for  calculating  fields  -  both  on-axis  and 
everywhere  -  along  with  numerical  ray  tracing  are  discussed  below  in  the  section  on  Numerical  Techniques. 

Another  issue  in  charged  particle  optics  are  those  mutually  repulsive  forces  that  charged  particles  exert 
between  themselves.  This  effect  is  known  as  space  charge  and  it  has  deleterious  effects  on  beam  quality. 
Traditionally,  the  beam  perveance,  IV~yi,  was  a  parameter  used  to  indicate  how  seriously  space  charge  will  affect  the 
beam  quality.  If  the  perveance  was  sufficiently  low  then  the  space  charge  effect  was  thought  to  be  negligible.  This 
was  the  result  from  work  done  on  vacuum  tubes  where  the  space  charge  was  modeled  as  a  continuum  of  charge 
density:  this  model  predicts  that  space  charge  acts  as  a  diverging  lens.  In  most  electron  beam  lithography  equipment 
the  values  of  perveance  are  such  that  this  lens  effect  is  negligible  and  so  space  charge  was  ignored  until  about  1970 
in  this  application.  However,  it  has  since  been  established  that  in  such  equipment  the  occurence  of  random 
interactions  is  sufficiently  frequent  to  cause  surprisingly  large  energy  spread  (the  Boersch  effect')  and  appreciable 
trajectory  bending  [3].  Because  of  the  random  nature  of  these  interactions  this  phenomenon  is  best  treated  using  a 
Monte  Carlo  technique  in  which  the  position  and  momentum  of  each  of  a  finite  population  of  particles  is  tracked  as 
the  panicles  move  down  the  column  [4],  A  whole  section  below  is  devoted  to  modeling  the  space  charge  effect 

Electron  sources  for  electron  beam  lithography  equipment  traditionally  have  been  relatively  simple 
structures  [5J.  However  point  sources  for  both  ions  and  electrons  pose  challenges  for  modelling  because  of  the  huge 
dispanty  between  the  nm-scale  geometries  of  the  tips  and  the  cm-scale  geometries  of  the  surrounding  electrodes. 

How  this  is  presently  tackled  is  described  in  the  section  on  Source  Modeling,  but  it  is  not  clear  that  a  satisfactory 
approach  yet  exists. 


Finally,  die  interaction  of  the  beam  with  the  solid  target  can  be  modeled  as  a  continuous  slowing  down  of 
the  particles  that  suffer  scanning  at  random  intervals  at  a  mean  frequency  set  by  scattering  cross  sections  of  the 
incident  particles  in  the  target  material.  Such  modelling  has  been  the  subject  of  countless  papers  since  1963  and  a 
brief  description  along  with  references  to  recent  literature  is  given  in  the  section  on  Beam-Target  interactions. : 


Numerical  Techniques 


With  the  exception  of  source  modeling  and  space  charge  effects,  much  of  the  computing  required  for 
Charged  Particle  sytems  goes  into  determining  electric  and/or  magnetic  fields  from  boundary  conditions  such  as 
potentials  on  electrodes  or  current  through  coils.  This  section  will  briefly  describe  some  of  the  numerical  techniques 
available  for  calculating  these  fields,  and  will  conclude  with  a  discussion  of  numerical  ray  tracing. 

For  a  cylindrical^  symmetric  system  (no  $  dependence),  any  potential  (electric  or  scalar  magnetic)  will 
have  the  following  Taylor  expansion  since  it  satisfies  the  Laplace  Equation  [6]: 

<D(r,  z)  =  <J>(z)  -  i<D"(z)rJ  +  -^-<I>""(z)r4  -  ...  (1) 

4  64 

This  says  that  the  potential  at  any  point  (r,  z)  can  be  determined  merely  by  knowing  <I>(r  =  0,  z) .  So,  in 

principle,  you  do  not  need  to  explicitly  compute  the  potentials  everywhere;  you  only  need  to  compute  them  along 
the  z-axis. 

Although  the  above  Taylor  expansion  will  only  be  valid  for  small  values  of  r,  much  useful  information 
may  still  be  extracted.  For  example,  the  paraxial  ray  equation  and  third  order  aberration  theory  can  be  derived  from  it. 
See  the  references  in  [7]  for  more  details.  In  addition,  a  modified  form  of  equation  (1 )  is  used  when  designing 
deflectors,  as  will  be  discussed  below. 

It  also  may  happen  that  one  needs  the  full  solution  to  the  potential/field  everywhere  inside  the  system,  not 
just  near  the  optical  axis.  This  will  be  required,  for  instance,  to  do  full  ray  tracing  of  particles  down  the  column  and 
especially  to  do  ray  tracing  near  an  electron  gun.  In  this  case,  the  full  partial  differential  equation  (Laplace’s 
equation)  must  be  solved. 


Discussed  below  are  the  three  most  popular  techniques  for  solving  for  potentials,  either  oo-axis  or 
everywhere  inside  the  system.  These  methods  «re  die  boundary  element  method  (BEM),  the  finite  difference  method 
(FDM),  sad  the  finite  element  method  (FEM). 


The  boundary  element  method,  also  known  as  the  charge-density  method,  is  one  way  to  quickly  determine 
$(r»  z)  fix’  electrostatic  problems.  What  one  does  is  to  break  up  the  electrodes  into  smaller  elements  (e.g.  an 

might  be  broken  up  into  sections  of  rings)  and  then  determine  what  charge  density  must  be  uniformly 
distributed  throughout  each  element  in  order  to  give  rise  to  the  prescribed  potentials  on  all  boundary  elements.  Here 
is  an  example: 


Annulus  at  5  V 


Annulus  at  0  V 


This  determination  of  the  charge  density  will  involve  solving  a  matrix  equation.  The  potential  at  any  point  in  space 
due  to  a  single  element  is  then  simply  found  by  using  the  formula 


Note  that  the  charge  density  p  is  uniform,  so  it  may  pulled  out  of  the  volume  integral;  typically  the  shapes  used  for 
the  boundary  elements  are  such  that  analytic  expressions  exist  for  the  remaining  volume  integral.  The  total  potential 
at  any  point  in  space  is  then  the  sum  of  the  potentials  from  all  elements.  Consult  the  references  in  (8]  for  more 
details. 


The  finite  difference  method  is  very  simple  to  understand.  Recall  that  the  nKopreudoo  of  the  first 
derivative  of  a  curve  is  that  his  the  slope  of  the  tangent  line  at  that  point  If  a  curve  c(x)  is  sampled  at  a  finite 
number  of  points,  say  a  grid  of  points  equally  spaced  by  the  distance  Ax,  then  we  might  approximate  the  fint 
derivative  of  c(x)  by 

dc(x).  c(a  +  Ax)  -  c(a)  f3) 

dx  |-  “  Ax 

Similarly,  the  second  derivative  may  be  approximated  by 

d2c(x),  c(a  +  Ax)  -  2c(a)  4  c(a  -  Ax)  //1X 

dx2  U  Ax2 

Similar  formulas  hold  for  partial  derivatives.  Then  any  second  order  PDE  with  boundary  conditions  should  be 
solvable  as  follows:  lay  out  a  mesh  of  grid  points  with  the  appropriate  dimensionality  (say  a  rectangular  grid  aligned 
with  the  x  and  y  axis  for  a  2D  problem)  and  insist  that  the  PDE  (as  expressed  with  the  above  approximations)  be 
satisfied  at  every  interior  mesh  point,  and  the  boundary  conditions  be  satisfied  at  every  boundary  mesh  point.  This 
will  generate  a  set  of  simultaneous,  linear,  algebraic  equations  (i.e.  a  matrix  equation)  which  can  be  solved  yielding 
the  potential  everywhere. 

The  finite  element  method  is  less  intuitive.  It  turns  out  that  solutions  of  PDEs  are  also  solutions  of  an 
associated  variational  problem.  What  is  important  is  that  finite  element  methods  are  similar  to  finite  difference 
methods  in  that  one  must  specify  a  grid  of  points.  One  also  must  specify  how  the  solution  is  assumed  to  vary  over 
a  given  mesh  surface  element;  a  first  order  FEM  assumes  the  solution  is  piecewise  linear  while  a  second  order 
solution  assumes  the  solution  is  piecewise  quadratic.  Finite  element  methods  are  more  general  than  finite  difference 
methods  (indeed.  FDM  can  be  shown  to  be  a  subcase  of  FEM).  One  consequence  of  this  is  that  FEM  can  be  easily 
applied  to  nonrectangular  meshes;  triangular  meshes  are  very  common.  This  may  be  a  large  advantage  in  accuracy 
when  the  boundary  is  very  nonrectangular  -  triangular  meshes  can  easily  be  chosen  to  have  the  boundary  mesh  sides 
lie  closely  along  the  boundary  [9].  Rectanguier  meshes  suffer  in  that  their  boundary  mesh  sides  often  weave  through 
a  rough  bondary  instead  of  lie  along  it.  In  addition,  FEM  is  the  only  way  to  handle  saturated  magnetic  lenses. 
Standard  algorithms  like  Gaussian  elimination  may  be  completely  inadequate  for  solving  the  huge  matrix  equation 
generated  by  FEM  (e.g.  solving  Poisson's  equation  using  second  order  FEM) .  Workers  in  the  field  have  now  begun 
to  employ  die  incomplete  Choleslu  conjugate  gradient  (1CCG)  method  [10]  and  report  considerable  success  with  it  - 
see  [1 1]  and  [12]. 


For  more  deoils  concerning  FDM  rad  FEM,  a  modem  reference  is  Hall  rad  Porschingf  1 3];  consult  pp.  1 59- 
183  of  Hawkes  [9]  for  die  application  of  these  methods  to  charged  particle  optics.  Other  references  to  the  application 
of  FEM  for  solving  electrostatic  fields  may  be  found  in  [14]. 

To  compare  BEM  with  FEM  for  solving  an  electrostatic  problem,  note  that  with  BEM  you  only  break  up 
the  electrodes  into  small  pieces  whereas  with  FEM  the  whole  space  is  broken  into  grid  elements.  This  means  that 
the  size  of  die  matrix  equation  which  must  be  solved  for  a  BEM  solution  is  significantly  smaller  than  that  for  FEM, 
which  is  a  significant  advantage  for  the  BEM  method.  On  the  other  hand,  with  FEM,  once  you  have  solved  the 
matrix  equation  you  have  the  solution  for  the  potential  everywhere  whereas  BEM  only  yields  the  charges  on  the 
electrodes.  To  get  the  potential  at  a  given  point  from  a  BEM  solution  you  must  add  up  the  potentials  due  to  the 
charge  on  each  boundary  element;  this  may  be  a  modest  computation  as  the  geometrical  integral  from  equation  (2) 
may  be.  say,  an  elliptic  integral  if  the  boundary  element  is  a  ring  (Reneau,  [8]).  So,  if  there  are  many  boundary 
elements  or  if  you  are  ray  tracing  a  sufficiently  large  number  of  particles,  the  FEM  may  be  more  economical.  One 
exception  may  be  ray  tracing  using  the  paraxial  ray  equation.  In  this  case,  since  you  only  need  the  potential  on-axis, 
the  BEM  may  be  more  economical. 


Given  the  solution  for  the  potential,  one  must  still  extract  the  fields  if  one  desires  to  perform  numerical  ray 
tracing.  Since  the  field  is  the  gradient  of  the  potential,  this  means  that  you  must  have  the  spatial  partial  derivatives 
of  the  potential.  If  a  grid  technique  like  FEM  has  been  used,  die  potential  is  only  known  on  the  grid  nodes  which 
means  that  an  accurate  method  for  smooth  interpolation  between  nodes  must  be  developed.  Chapter  13  of  Hawkes 
[9]  goes  into  detail  describing  some  of  the  available  methods;  see  also  some  of  the  references  in  [14],  Lunney  et  al 
[15]  have  developed  a  new  method  based  on  multipole  expansion  which  they  claim  is  a  superior  alternative  to 
interpolation  based  schemes  for  computing  the  field. 


Once  the  fields  have  been  determined,  an  accurate  method  for  ray  tracing  the  particles  must  be  employed.  If 
a  cylindrically  symmetric  system  is  under  consideration  and  one  has  solved  for  the  axial  field,  then  one  may  elect  to 
solve  for  the  trajectories  using  the  paraxial  ray  equation  [7].  For  a  magnetic  lens  system,  the  paraxial  ray  equation  is 

-  o 


where  q  is  the  charge  to  mass  ratio  and  Vr  is  the  reiativisticaliy  corrected  beam  voltage.  It  is  just  an  ordinary 
differential  equation  for  the  radial  coordinate  of  the  trajectory  as  a  function  of  the  axial  coordinate,  so  the  standard 
methods  such  as  Runge-Kutta  or  Predictor-Corrector  may  be  employed  to  numerically  solve  it.  See  Press  et  al  for  a 
good  description  of  these  techniques  [16]. 


If  the  accuracy  of  the  paraxial  ny  equation  is  insufficient,  tben  more  precise  techniques  must  be  employed. 
One  such  accurate  way  of  rayjnciag  is  to  consider  a  Taylor  aeries  expansion  io  oitler  to  determine  tbe  future  position 
of  a  particle  given  its  present  position: 

r(t  +  At)  *  r(t)  +  r'(t)Ai  +  ^r"(t)At2  +  ^r~(t)Al3  +  ...  (6) 

2  6 

where  primes  denote  differentation  with  respect  to  the  time  t  If  v(t)  is  defined  as  the  velocity  and  a(t)  as  the 
acceleration,  we  can  rewrite  equation  (6)  as: 

r(t  +  At)  a  Tit)  +  v(t)At  +  |a(t)At2  +  |a'(t)Ai3  +  ...  (7) 

2  6 

Making  use  of  Newton's  Second  law  (F  =  ma),  we  can  express  this  as: 

r(t  +  At)  =  r(t)  +  v(t)At  +  -^-F(r(t),  t)At2 

2m 


~F(v(t)«  V)F(r(t),  t)  +  v1  t)~j At3 
)m  L  ot 


The  force  F  for  electromagnetic  forces  is 

F(r(t),  t)  =  q(E  +  vxB)  (9) 

and  is  in  general  a  function  both  of  particle  position  and  time  since  E,  v,  and  B  in  general  depend  on  position  and 
time.  From  equation  (8)  you  can  trace  the  particles'  trajectory  by  advancing  the  time  in  small  increments  At.  This 
is  known  as  the  power  series  method  of  ray  tracing.  The  terms  shown  in  equation  (8)  constitute  a  third  order  power 
series  method;  it  is  common  to  drop  the  last  term  and  do  a  second  order  method  if  less  accuracy  will  suffice. 

An  alternative  to  the  power  series  approach  would  be  to  numerically  solve  the  differential  equation  for  the 
trajectory.  This  differential  equation  is  just  Newton's  Second  Law: 

~ =  — ~F(r(t),  t)  (10) 

at  m 

where  F  is  again  given  by  equation  (9).  It  is  simple  to  see  that  in  rectangular  coordinates.  Otis  equation  breaks 
down  into  three  coupled  ordinary  differential  equations.  Then  the  standard  techniques  like  Runge-Kutta  or  Predictor- 
Corrector  may  be  used  to  solve  them  ( 16];  see  [  17}  as  an  example  in  the  literature  where  this  is  applied  to  charged 
panicle  trajectories. 


A  very  accurate  technique  for  tracking  a  panicle's  motion  in  a  force  field  is  the  technique  of  Nystrom 
integrators,  developed  in  1925  by  Nystrom  {18].  This  technique  is  also  very  complicated  which  has  hindered  its 
acceptance.  Lear  has  exploited  the  capabilities  of  computers  to  ease  the  difficulty  of  implementing  the  method  [19]. 
His  application  was  actually  the  study  of  orbital  motion,  but  perhaps  there  may  be  application  in  charged  panicle 
optics. 


Source  Modeling 


There  are  several  issues  involved  in  modeling  charged  particle  sources.  First,  there  must  exist  a  good 
physical  model  for  the  emission  process.  For  thermionic  electron  sources,  the  Richardson-Dushman  equation  with 
Schottky's  field  enhancement  correction  has  proven  to  be  a  good  modeL  Held  emission  electron  sources  are 
adequately  modeled  by  the  Fowler-Nordbeim  equation.  A  difficult  source  to  model  is  liquid  metal  ion  sources, 
because  the  tip  is  not  fixed  but  changes  shape  during  the  emission  process  in  a  m»wnw  which  is  hard  to  predict 

Once  the  appropriate  model  has  been  selected,  a  source  modeling  program  needs  to  very  accurately  determine 
the  fields  near  the  source.  This  is  both  because  the  emission  process  itself  always  strongly  depends  on  the  field 
configuration  and  also  because  the  electrons  will  be  ray  traced  upon  emission.  The  most  common  candidate  for 
determining  the  fields  are  die  grid  methods  (FDM  or  FEM).  In  a  typical  electron  source,  the  emission  area  may  be 
as  small  as  a  few  hundred  angstoms  while  the  extraction  electrodes  may  be  spaced  millimeters  away  from  the  source. 
This  huge  length  scale  difference  means  that  a  grid  method  for  solving  the  fields  cannot  have  constant  grid  spacing  - 
this  would  require  too  many  elements  to  adequately  model  the  source  region.  Thus,  the  grid  must  change  scale  from 
small  grid  spacing  near  the  source  to  large  grid  spacing  near  the  eiccdodes.  This  is  almost  always  far  too  tedious  for 
a  human  to  specify  the  grid  sizing  by  hand,  so  a  good  source  modeling  program  should  include  a  routine  for 
automatically  generating  a  grid  with  adaptive  length  scales. 

Lastly,  it  is  usually  vital  to  include  space  charge  effects  when  modeling  the  emission  process  -  the 
repulsion  caused  by  electrons  which  have  already  been  emitted  can  greatly  affect  the  emission  of  additional  electrons. 
Furthermore,  the  presence  of  other  electrons  will  modify  the  trajectories  that  would  result  if  only  fields  from  the 
electrodes  were  present.  See  below  for  a  whole  section  which  is  devoted  to  techniques  for  modeling  space  charge 
effects. 


One  of  the  most  comprehensive  solutions  to  modeling  electron  sources  is  the  program  SOGUN  reported  by 
Zhu  and  Munro  [12].  Their  program  can  handle  either  thermionic  or  field  emission  guns.  They  use  a  2nd  order, 
isoparametric  FEM  code  to  solve  for  the  electric  potential.  The  sides  of  the  grid  elements  in  this  technique  arc  not 
lines  but  are  quadratic  curves  -  this  enables  the  grid  lines  to  almost  perfectly  conform  to  the  boundary  shapes  (e  g. 
curved  cathodes  and  electrodes).  The  fact  that  this  is  a  2nd  order  method  means  that  the  potential  can  be  solved  with 
higher  accuracy.  Once  the  potential  has  been  solved  for,  the  electric  field  needs  to  be  extracted.  This  is  one  difficulty 
with  isoparametric  grids,  but  Zhu  and  Munro  report  a  new  algorithm  to  do  this.  The  commercial  version  of  this 
program,  SOURCE,  also  includes  automesh  generation  so  that  logarithmic  changes  in  grid  size  -  which  is  needed  for 
field  emission  guns  -  can  be  done  by  computer  [20].  Ray  tracing  is  done  by  a  third  order  power  series  method  [21], 


and  apace  charge  effect*  are  modeled  by  iteratively  aolvinf  Poisson's  equation  as  beiow.  The  program  is 

abk  to  determine  the  total  boon  current,  the  crossover  poaitioa  and  size,  and  beam  aberrations.  Sample  outputs 
from  Mumo's  software  may  be  found  on  the  next  two  pages,  where  they  model  a  thermionic  LaB6  gun  and  a  fi»M 

emission  gun. 

Other  earlier  work  on  modeling  sources  may  be  found  in  Weber  [22],  Hemnannsfddt  [23],  K»«»e  [24],  and 
Renau  [25].  Browning  [26]  has  successfully  employed  the  BEM.  See  Hawkes  [27]  for  a  discussion  of  an  analytical 
model  for  source  region  space  charge  effects.  A  reference  to  modeling  liquid  metal  ion  sources  may  be  found  in  Cui 
and  Tong  [28].  Mohammed  and  Garcia  [29]  discuss  automesh  generation  for  electrostatic  problems. 


Lens  Modeling 


Lens  modeling  is  (be  most  well  developed  area  for  computer  modeling  of  charged  particle  systems,  with 
early  programs  appearing  over  20  yen  ago  [2]. 

The  common  feature  of  almost  all  probe  fbnninf  leases,  whether  electrostatic  or  magnetostatic,  is  that  they 
have  cylindrical  symmetry.  For  example,  a  typical  electrostatic  tom  comitt  of  ring  electrodes  and  a  typical 
magnetostatic  lens  is  cylindrical  current  coils  surrounded  by  high  permeability  materials  which  concentrate  the 
magnetic  field  across  a  small  gap.  This  symmetry  reduces  a  3D  field  problem  to  a  2D  problem. 

The  historical  method  of  handling  this  problem  is  to  consider  a  Taylor  series  exapansion  of  the  potential 
(Equation  1).  This  method  requires  knowledge  of  the  potentials  on-axis,  which  can  be  obtained  by  any  of  the 
described  in  the  Numerical  Techniques  section.  Once  these  have  been  obtained,  the  on-axis  fields  (electric  or 
magnetic)  are  easily  determined  and  then  the  paraxial  ray  equation  -  an  ordinary  differential  «qn«t«n«i  (see  5) . 

must  be  numerically  solved  for  certain  trajectories  (the  "principal  rays").  These  trajectories  are  then  used  in 
numerical  integrations  to  obtain  the  aberration  coefficients;  important  ones  are  spherical  and  chromatic  aberrations. 
These  aberration  coefficients  are  then  used  to  predict  lens  performance.  Chu  and  Munro  [7]  report  on  a  program 
capable  of  modeling  both  lenses  and  deflectors  in  this  manner. 

We  give  a  worked  example  below  of  how  to  model  the  performance  of  a  magnetic  lens.  This  will  give  the 
reader  an  idea  of  what  computer  modeling  of  charged  particle  optics  involves.  This  particular  example  is  taken  from 
a  manual  Munro  has  provided  our  lab  [30]. 


The  leu  geometry  we  will  consider  is  shown  on  the  next  pegs,  labeled  as  Tig.  18b”.  The 
first  task  win  be  to  determine  the  axial  magnetic  field  (also  known  as  the  flux  density),  using 
Monro's  program  Mil.  Because  of  die  symmetry  of  die  leu,  we  need  only  solve  over  the 
region  AYZD;  this  region  is  shown  in  detail  in  "Fig.  19b*.  Mil  expects  that  you  have 
divided  the  geometry  into  quadrilaterals,  as  shown  in  "Fig.  20b”,  in  which  edges  of  the  pole 
piece  lie  along  edges  of  die  quadrilaterals.  You  will  need  to  teU  Ml  1  the  exact  location  of 
these  quadrilaterals.  Further  more,  since  Ml  1  uses  a  FEM  technique  to  solve  for  the  magnetic 
scalar  potential  everywhere  inside  AYZD,  you  will  need  to  tell  Ml  1  what  density  of  mesh 
lines  to  use  inside  each  quadrilateral.  This  is  sketched  out  in  Tig.  21b”:  the  numbers  along 
the  lines  YZ  and  ZD  are  the  locations,  in  millimeters,  of  the  corresponding  quadrilateral 
vertex:  the  numbers  along  the  lines  AD  and  AY  label  the  mesh  lines  thru  the  quadrilateral 
edges.  The  completed  mesh  is  shown  in  "Fig.  22b”.  Now  the  steps  outlined  above  are 
something  that  you  die  user  do  on  paper  -  the  actual  input  data  supplied  to  the  computer,  in 
the  form  of  a  text  file,  is  shown  in  "Table  lb”.  This  is  the  format  of  "Table  lb”:  the  first  two 
big  blocks  of  numbers  are  die  axial  ("z")  and  radial  ("r")  coordinates  of  the  quadrilateral 
vertices;  the  numbers  1-7-12-22-27  and  1-5-1B-15-20  on  the  top  and  left  of  these  blocks  are 
the  mesh  line  labels  (compare  with  "Fig.  21b").  The  next  line  of  numbers,  1-22-5-15-1000 
tells  the  location  of  the  pole  piece  and  its  relative  permeability.  The  final  two  columns  of 
numbers  tell  the  boundary  potential  distribution  (in  amptums).  The  output  of  this  program  is 
the  axial  magnetic  field,  shown  in  "Table  2b". 

Now  that  the  axial  magnetic  field  has  been  obtained,  we  may  go  on  to  examine  the  objective 
properties  of  the  lens.  Munro's  program  M21  will  calculate  the  excitation  parameter  and  - 
depending  on  the  magnification  condition  •  the  object  and/or  image  plane,  the  objective 
principal  plane,  the  objective  focal  length,  the  objective  magnification,  the  spherical 
aberration  coefficient,  the  chromatic  aberration  coefficient,  and  the  magnetic  field  at  the  object 
and/or  image  plane.  These  are  the  main  properties  needed  to  predict  lens  performance.  Along 
with  the  data  file  containing  the  axial  magnetic  field,  the  user  must  also  input  to  M2I  the 
initial  beam  voltage,  the  increment  in  beam  voltage,  the  number  of  beam  voltages,  the 
magnification  condition  (zero,  low,  high,  or  infinite),  and  possibly  the  position  of  the  object 
or  image  plane  (if  the  magnication  is  either  low  or  high).  Examples  of  the  output  of  M2 1  for 
the  lens  of  Tig.  18b"  are  found  in  "Table  14"  and  "Table  15". 
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0.0001 

15000. 

4.08  -200.00 

35.67 

0.179 

685.67 

39.69 

0.0001 

16000. 

3.95  -200.00 

38.35 

0.192 

856.09 

43.28 

0:0001 

17000. 

3.83  -200.00 

41.08 

0.206 

1059.17 

47.03 

0.0000 

18000. 

3.73  -200.00 

43.88 

0.220 

1299.78 

50.94 

o.oooo- 

19000. 

3.63  -200.00 

46.75 

0.234 

1583.48 

55.03 

0.0000 

20000. 

3.54  -200.00 

49.68 

0.249 

1916.48 

59.29 

0.0000 

THE  ABOVE  RESULTS  ARE  FOR-  AN  EXCITATION 
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As  noted  in  the  Introduction,  the  above  technique  of  Taylor  series  expansion  only  handle  third  order  effects 
and  this  may  result  in  poor  accuracy  for  many  situations.  While  the  theory  can  be  extended  to  higher  orders,  the 
resulting  equations  are  much  more  complicated  so  that  an  alternative  technique  is  necessary.  One  such  alternative 
method  of  lens  modeling  is  to  give  the  lens  design  parameters  (e.g.  electrode  shapes  and  potentials)  to  a  program 
which  determines  the  fields  everywhere  in  the  lens.  Direct  ray  tracing  can  then  be  used  to  see  bow  the  electrons 
navel.  This  easily  yields  characteristics  such  as  optimal  focal  plane  position  and  focal  spot  size  with  high  accuracy. 
The  fields  are  determined  by  first  solving  for  the  potential,  typically  by  either  the  FDM  or  FEM  mesh  methods.  A 
recent  innovation  is  the  use  of  the  ICCG  method  for  solving  the  associated  matrix  equation  for  FEM  determination 
of  the  potential;  see  Lencova  and  Lenc  [1 1].  Lencova  and  Wisaeiink  [31]  also  introduced  automeshing  capabilites  for 
their  lens  design  program  -  although  there  is  a  less  serious  need  for  this  capability  in  lens  design  as  compared  to 
source  modeling.  Any  good  program  for  designing  magnetic  lenses  should  have  the  capability  of  handling  magnetic 
saturation  effects. 

A  commercial  software  package  which  is  capable  of  modeling  charged  particle  lenses  using  this  second 
technique  of  direct  ray  tracing  is  the  program  OPTICS  put  out  by  Munro's  company  [20].  A  typical  example  of  the 
capability  of  Munro's  software  is  shown  in  the  top  picture,  (a),  on  the  next  page.  The  near  vertical  lines  are  the 
equipotentials  between  the  two  electrodes.  The  near  horizontal  lines  are  trajectories  for  electrons  with  different 
starting  heigths.  The  focal  plane  and  spot  size  are  easily  determined  from  where  the  trajectories  are  most  closely 
crossed  together.  A  plot  of  focal  position  vs.  square  of  ray  slope  is  on  the  bottom  -  note  the  fifth  order  effect  which 
the  above  third  order  theory  would  miss. 


Examples  of  earlier  work  in  computer  modeling  of  electrostatic  lenses  may  be  found  in  the  program 
CIELAS  by  Hill  and  Smith  [32]  and  the  program  ELOP-GELOP  by  van  Oostnim  [33]. 


Deflector  Modeling 


Deflector  is  timiUr  to  leas  modeling  in  that  one  mutts  to  determine  the  field  distribution  caused 

by  electrode  or  current  coil  arrangements  and  then  use  electron  ray  tracing  or  aberration  coefficients  to  predict  the 
performance.  The  main  difference  between  deflectors  and  lenses  is  that  deflectors  never  have  cylindrical  symmetry,  so 
determining  the  field  distribution  is  a  3D  problem.  In  principal,  one  could  use,  say,  a  3D  FEM  program  to  solve  for 
the  fields.  In  practice,  3D  models  are  so  much  more  demanding  of  computer  time  and  memory  capacity  compared  to 
2D  models  that  they  have  not  yet  been  used  for  deflector  design.  Instead,  approximate  methods  have  been  developed 
to  cope  with  the  3D  problem. 

It  turns  out  that  electrostatic  deflectors,  while  lacking  cylindrical  symmetry,  typically  still  possess  some 

degree  of  axial  rotation  symmetry.  In  particular,  the  most  common  electrostatic  deflectors  have  invariance  with 

respect  to  rotations  of  either  90  or  45  degrees.  In  such  a  case,  the  3D  field  distribution  can  be  calculated  via  Fourier 

expansion  as  a  sum  of  harmonic  terms: 

0(r,  0,  z)  =  <J>,(r,z)cos(0)  +  <J>3(r,  z)cos(30)  +  ... 

=  £-fj(z)r  +  if','(z)r3 

(the  rotational  symmetry  eliminates  all  even  harmonic  terms).  The  common  approximation  is  to  neglect  the  higher 
order  terms  and  only  deal  with  terms  shown  in  the  equation  above;  this  may  be  an  inadequate  approach  for  large  angle 
deflectors.  Methods  for  calculating  either  the  functions  or  the  functions  f]  and  f3  for  both  magnetic  and 

electrostatic  deflectors  are  outlined  in  Munro  and  Chu  [34].  These  methods  include  invoking  the  Biot-Savart  law  for 
magnetic  deflectors  if  no  ferromagnetic  materials  are  present  -  otherwise  a  FEM  code  is  run.  Likewise,  electrostatic 
deflectors  can  be  solved  either  by  the  BEM  method  or  FEM.  Lencova  [35]  reports  a  code  for  handling  tapered 
magnetic  deflectors. 

There  are  many  lenses,  such  as  multipole  lenses,  slit  lenses,  grid  lenses,  and  concentric  hemispherical 
analysers  which  are  not  cylindrically  symmetric.  Multipole  lenses,  which  are  known  to  be  capable  of  correcting 
beam  aberrations  [36],  can  be  modeled  exactly  like  deflectors  if  they  have  the  same  symmetry  properties;  Smith  and 
Munro  report  a  multipole/deflector  program  [37],  The  other  types  of  lenses  may  be  modeled  in  a  manner  similar  to 
equation  (11).  except  that  the  lack  of  symmetry  will  mean  that  the  even  harmonic  terms  will  need  to  be  included. 

Future  trends  for  deflector  modeling  will  probably  be  to  increase  the  accuracy  of  the  field  calculations  and  to 
model  traveling  wave  deflectors. 
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Tolerances 


A  vital  tool  for  charged  particle  optics  design  is  a  program  which  can  «i*nt tolerances  -  the  maximum 
amount  of  mechanical  imperfections  which  do  not  move  the  performance  oat  of  specification.  Typical  imperfections 
are  things  like  errors  in  machining  or  lens  elements  not  being  property  aligned. 

The  basis  for  all  tolerance  modeling  programs  is  Stunock's  Principle  [38].  Sturrock's  Principle  says  that 
the  effect  of  moving  a  point  P  on  an  electrode  by  the  vector  Sfp  is  equivalent  to  leaving  the  point  unchanged  but 

changing  the  potential  by  84>p  =  -  V<I>p  •  8rp.  This  is  really  convenient'  to  model  the  effect  of,  say,  a  lens  . 

being  machined  elliptical  instead  of  cylindrical,  one  could  actually  model  the  lens  as  being  cylindrical  but  now  have  a 
changed  boundary  potential  (via  Stunock's  principle).  The  fact  that  the  boundary  potential  is  no  longer  cylindrically 
symmetric  means  that  the  same  techniques  developed  for  deflector  design  must  now  be  invoked,  namely  Fourier 
expansion  as  a  sum  of  harmonic  terms.  Note  that  unlike  deflector  design,  the  even  harmonic  terms  may  be  nonzero. 
This  will  now  give  you  an  approximation  for  the  axial  potential  C>(r  =  0,  z)  which,  when  combined  with  a 

unified  lens/deflector  aberration  theory  [7],  allows  you  to  compute  the  aberrations  by  numerically  solving  integrals. 

Liu  [39]  and  the  references  they  cite  describe  a  modem  approach  to  lens  and  deflector  toierancing  software. 
Arc  hard  [40]  is  a  useful  reference  to  assist  in  understanding  Sturrock's  original  paper. 


Space  Charge  Modeling 


All  charged  particle  optical  systems  suffer  from  space  charge  effects.  This  refers  to  the  mutually  repulsive 
electric  forces  that  similarly  charged  panicles  exert  on  each  other  [41].  Hence,  a  beam  of  electrons  will  actually  have 
different  trajectories  when  going  through  a  lens  than  those  that  would  be  predicted  from  ray  tracing  when  only  the 
static  fields  from  the  lens  are  taken  into  account.  These  mutually  repulsive  forces  have  deleterious  effects;  the  axial 
components  cause  the  beam  energy  spread  to  widen  (the  "Boersch  effect")  and  the  radial  components  cause 
broadening;  both  degrade  the  focusing  ability  of  the  system. 


With  an  estimate  for  the  beam  current  density,  a  first  order  model  for  the  effects  of  space  charge  is  to 
assume  that  die  charge  density  is  static  and  continuosly  distributed  over  the  column  and  then  to  solve  for  the  new 
fields  in  the  system.  For  instance,  to  model  space  charge  effects  in  an  electrostatic  lens  in  this  manner,  one  would 
be  solving  Poisson's  equation  instead  of  the  simpler  Laplace's  equation  (which  is  used  when  only  electrodes  are 
considered).  Fortunately,  this  only  requires  slight  modification  of  the  codes  used  in  lens  design. 

One  way  to  obtain  an  estimate  for  the  current  density  is  to  calculate  the  fields  with  no  charge  assumed 
present,  trace  the  trajectories  of  a  bunch  of  particles  through  this  system  (where  the  particles  start  out  with  random 
positions,  energies,  and  entry  times  -  the  randomness  is  a  distribution  obtained  from  some  model),  and  then  to  use 
these  trajectories  to  estimate  the  charge  density.  Now  one  can  solve  Poisson's  equation  for  a  more  accurate  solution. 
One  could  continue  to  iterate  on  the  charge  densities  derived  from  successive  solutions  if  desired.  This  is  die  method 
used  by  Munro  in  his  electron  gun  program  to  simulate  space  charge  effects  near  the  source.  He  has  found  that  this 
method  converges  to  a  solution  in  about  10  iterations  (see  p.  1866  of  reference  [12]). 

A  more  effective  method  to  model  space  charge  is  by  full  Monte  Carlo  simulation.  In  such  a  simulation, 
particles  again  are  assigned  random  starting  positions,  energies,  and  entry  times  according  to  some  distribution 
model.  Then  direct  ray  tracing  is  performed,  using  both  die  static  fields  from  the  optical  elements  and  the  dynamic 
electric  fields  generated  by  all  the  pair  interactions.  This  method  is  appealing  because  it  correlates  with  what  is 
happening  physically  as  a  particle  beam  moves  down  a  column.  There  are  effects  predicted  by  these  Monte  Carlo 
studies  which  the  above  first  order  models  do  not  account  for.  In  particular,  when  two  real  electrons  m»ire  a  close 
approach  they  experience  higher  angle  scattering  than  the  scattering  predicted  from  a  continuum  distribuuon.  This 
effect  is  known  as  residual  stochastic  Coulomb  interactions.  See  Hawkes'  book  [42]  for  a  brief  discussion 

Numerically,  there  are  two  issues  to  be  faced  when  implementing  a  Monte  Carlo  simulation.  First,  if  there 
are  N  charged  particles  in  the  column,  then  the  number  of  pair  interactions  which  must  be  ralrnlatwi  grows  as  N2. 
This  makes  a  full  Monte  Carlo  simulation  of  space  charge  effects  computationally  expensive  when  many  particles 
are  present  (e.g.  at  high  current  levels).  Unfortunately,  this  is  precisely  when  you  are  most  concerned  about  space 
charge  effects.  Second,  the  accuracy  of  the  numerical  ray  tracing  method  is  highly  dependant  upon  choosing  small 
time  steps  when  two  particles  move  close  to  each  other.  So,  an  appropriate  algorithm  must  be  employed  to 
dynamically  adjust  the  time  step.  This  is  additional  overhead  and  the  smaller  time  steps  required  also  slow  the 
program  down.  Other  than  these  two  considerations,  a  space  charge  modeling  program  is  easy  to  implement  as  the 

electric  field  from  pairwise  interaction  of  electrons  is  trivial  to  calculate  (the  electrostatic  expression  is 
E  =  e  f /(4ne0r3)). 

In  spite  of  the  computational  expense  involved  in  a  full  Monte  Carlo  simulation,  many  people  have 
developed  codes  employing  it.  Munro's  program  COULOMB  uses  it  [43];  his  program  can  model  Gaussian  round- 
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accaieratiag  or  retarding  sections.  ShnnoysaM  [44]  has  apparently  dooe  a  full  Mooie  Cario  investigation  imp  space 
charge  effects  in  electron  sources.  Groves  [45]  and  Sasaki  [46]  am  among  others  who  report  full  Monte  Carlo 
simulation  codes. 

In  order  to  speed  up  space  cfaarfe  calculations  while  keeping  some  of  the  benefits  of  a  full  Monte  Carlo  code  . 
(e.g.  the  hard  scattering  effects),  there  have  been  several  attempts  at  approximate  Moote  Carlo  methods.  One 
approach  is  that  of  Alice  [47],  in  which  they  restrict  the  number  of  pair  interactions.  Earlier  work  by  Yau  [48]  had 
investigated  the  possibility  of  having  a  fixed  sphere  of  influence  about  each  electron  such  that  only  nearest  neighbor 
interactions  are  considered.  AUee  incorporated  a  more  sophisticated  approach  in  which  the  size  of  this  sphere 
(actually  a  cube)  is  dynamically  adjusted  so  that  on  average  only  a  specified  number  of  electrons  are  inside  it.  The 
dynamic  "sphere"  adjustment  is  cleverly  implemented  so  that  only  a  few  comparisons  per  electron  need  be  done. 

They  found  that  using  only  about  10  nearest  neighbors  gave  results  comparable  to  including  all  the  pair  interactions. 

As  a  result  of  considering  fewer  interactions,  the  calculation  time  will  grow  as  N14,  which  is  superior  to  the  N2 
„  growth  for  a  fell  Monte  Carlo  simulation. 

Another  possibility  is  the  "Fast"  Monte  Carlo  technique  developed  by  Jansen  [49],  [SO].  In  this  method, 
the  particles  are  given  random  initial  positions  and  velocities  just  as  in  a  full  Monte  Carlo  simulation  The 
difference  is  that  the  ray  tracing  is  not  done  numerically  but  analytically:  if  one  can  assume  that  space  charge 
interactions  only  cause  small  deviations  in  trajectory,  then  analytic  formulas  can  be  derived  to  calculate 
deviations.  These  formulas  may  still  require  numerical  methods  to  solve,  but  this  turns  out  to  be  much  quicker  than 
numerical  ray  tracing.  Jansen  claims  a  speedup  factor  of  10  to  100  times  for  his  method.  The  assumption  about 
small  deviations  is  not  valid  for  all  beam  conditions  -  in  particular,  it  becomes  false  for  large  beam  currents. 
Nevertheless,  for  small  to  medium  beam  currents,  his  method  may  provide  the  fastest  way  to  get  a  good  estimate  for 
space  charge  effects. 


Special  mention  should  be  made  here  about  space  charge  modeling  for  focused  ion  beam<FIB)  systems. 
FIB*s  are  more  prone  to  suffer  from  space  charge  effects  than  electron  systems  for  two  reasons.  First,  ions  are  many 
thousands  of  times  heavier  than  electrons,  so  for  a  given  energy  they  move  much  slower.  This  means  that  to 
achieve  a  given  beam  current,  the  density  of  ions  must  be  much  greater  than  the  density  of  electrons  in  a  beam  with 
the  same  current.  Second,  FIB'S  are  typically  operated  at  higher  current  levels  than  electron  beams.  Another  way  to 
recognize  the  susceptibility  of  FIB  systems  to  space  charge  is  the  fact  that  a  tightly  focused  charged  particlebeam  is 
not  in  equilibrium,  so  the  longer  time  of  flight  means  greater  relaxation  towards  equilibrium;  see  [42]  and  [50].  Two 
workers  who  report  on  Monte  Carlo  codes  for  FIB  systems  are  Nanim  [5 1  ]  and  Vijgen  [52].  Nanim's  code  was 
capable  of  implementing  a  dynamic  sphere  of  influence  to  limit  the  number  of  panicle  interactions  and  dynamic  time 
steps  while  Vijgen's  modeling  used  Jansen's  "Fast"  Monte  Carlo  programs. 


Beam-Target  Interactions 


lithography.  PtriMBplt,  to  quantify  images  or  undtwtsad  charging  effects  hi  i  Scaadg  Sectrou  Microscope 
(SEM)  will  require  detailed  knowledge  of  the  emission  of  secondary  dectrous  [55].  Modeling  beam-target 
interactions  may  be  very  important  for  lithographic  application*  where  different  parameters  such  as  beam  energy  may 
greatly  affect  proximity  exposure  or  sample  damage. 

An  introduction  to  the  physics  of  electron  scattering  and  diffusion  in  solids  can  be  found  in  Retmer  [56]. 
David  C.  Joy  has  published  an  excellent  introduction  to  Monte  Carlo  electron  emulations  [57].  The 

typical  method  makes  several  simplifying  assumptions.  Hist,  it  is  assumed  that  elastic  scattering  can  be  sepmaied 
from  inelastic  scattering;  elastic  scattering  is  assumed  responsjhle  for  changing  the  dhection  of  the  electron  trajectory 
while  inelastic  scattering  is  assumed  to  only  cause  energy  loss.  Next,  for  elastic  scattering,  the  solid  lattice  of  atoms 
is  viewed  as  a  collection  of  die  individual  atoms.  Finally,  die  inelastic  scattering  is  assumed  to  continuously  dnm 
an  electron  of  its  enetgy  at  a  rate  given  by  the  Bethe  formula  [58].  The  elastic  scattering  cross  section  used  is  the 
screened  Rutherford  formula  for  an  isolated  atom;  this  formula  predicts  a  mean  free  pnlb  winch  is  used  to  determine 
in  a  random  fashion  how  long  an  electron  will  travel  until  it  next  scatters.  At  its  next  witffag  site,  the  scatter 
angle  will  next  be  determined  and  the  whole  process  will  be  repeated  until  the  electron  runs  out  of  energy  due  to  the 
continuous  inelastic  scattering  taking  place  between  elastic  scattering  events.  Several  thousand  compute  trajectories 
will  yield  an  estimate  with  a  few  percent  error  of  such  characteristics  as  electron  penetration  and  energy 

deposition.  Joy  has  made  his  programs  freely  available  to  the  public  [59]. 

One  big  weakness  of  dis  approach  is  the  use  of  the  screened  Rutherford  elastic  cross  section  formula.  It  is 
known  that  this  formula  has  limited  validity  -  in  particular,  it  is  bad  for  elastic  scattering  off  of  heavy  at 

entrgiks  \bek>w,  say,  10  keV).  A  mote  accurate  method  of  obtaining  the  elastic  cross  tfctwa  is  p-riil  wave 
expansion  using  the  Mott  scattering  formula,  but  this  method  requires  an  extensive  calculation  (several  hours  of 
supercomputer  time)  each  time  you  want  to  use  ii|60].  Browning  1611  has  determined  an  empirical  expression 
which  interpolates  over  tabulated  Mott  cross  section  data  for  a  large  range  of  electron  energies  and  atomic  aims.  This 
now  allows  the  significantly  more  accurate  Mott  cross  sections  to  be  used  in  a  practical  manner  for  Monte  Carlo 
scattering  simulations. 


A  more  difficult  area  to  accurately  model  is  inelastic  scattering.  There  are  several  important  inelastic 
scattering  mechanisms  like  atomic  inner  shell  ionization  (X-ray  and  Auger  electron  generation;  see  [62]),  scattering 
from  both  valence  and  conduction  electrons,  and  collective  phenomena  such  as  plasmon  excitations.  These  effects 
may  depend  not  only  on  atomic  constituents,  but  also  on  crystal  structure.  There  is  a  lack  of  accurate  experimental 
data  in  this  field.  David  Joy  has  compiled  a  database  containing  all  the  published  data  he  could  find  on  electron 
backscatter  coefficients,  secondary  electron  yield,  stopping  power,  and  X-ray  ionization  cross-sections  [63].  A  glance 
at  the  data  reveals  rather  astonishing  discrepancies  between  the  results  reported  by  different  authors  and  large  gaps 
between  data  points  for  many  elements.  The  experiments  done  to  gather  this  data  are  particularly  difficult  because 
the  experimental  apparatus  is  effectively  put  of  the  sample;  no  one  has  yet  done  a  complete  analysis  of  how  to 
subtract  this  effect  out  [64].  As  a  result  of  this,  it  is  not  at  all  uncommon  to  see  researchers  publish  papers  in 
which,  say,  the  forward  and  back  scattered  currents  do  not  add  up  to  100%  of  the  incident  beam.  In  addition,  the 
experimentalists  have  been  irresponsible  by  not  fully  reporting  the  sample  characteristics.  For  example,  it  might 
make  a  big  difference  in  the  data  whether  a  carbon  sample  is  amorphous,  graphite,  or  diamond.  The  topography  of 
the  sample  also  has  a  large  effect  on,  say,  secondary  electron  yield  (trenches  will  act  as  Faraday  cups)  yet  this  is 
another  characteristic  that  may  go  unreported  [64]. 

References  to  earlier  work  cm  Monte  Carlo  studies  of  electron  wararing  in  solids  may  be  found  in  [65]; 
references  which  focus  on  electron  scattering  as  applied  to  electron  beam  lithography  may  be  found  in  [66]. 

There  is  also  interest  in  ion  scattering  in  solids.  People  have  done  Monte  Carlo  studies  of  ion  scattering  in 
solids,  but  the  problem  is  much  more  challenging  than  electron  scattering.  The  reasons  for  this  include  atom  knock¬ 
out  effects  (lattice  atoms  may  get  knocked  out  of  their  sites  and  in  turn  knock  out  other  atoms  causing  a 
effect),  crystal  effects  like  channeling  (ions  may  travel  further  along  certain  crystal  orientations),  and  the  fact  that  a 
heavy  enough  ion  dose  will  modify  the  characteristics  of  the  sample  (i.e.  turn  it  amorphous).  The  people  who  are 
concerned  with  modeling  ion  scattering  in  solids  appear  to  be  mainly  interested  in  ion  implantation  processes  for  the 
semiconductor  industry;  see  the  references  in  [67]  for  more  details. 


Future  directions 


As  this  paper  has  shown,  existing  programs  model  many  topics  in  charged  particle  systems.  However, 
almost  all  of  these  programs  exist  in  isolation  so  that  one  big  feature  would  be  the  concatenation  of  many  of  the 
existing  programs  into  one  easy  to  use  program.  For  instance,  while  there  are  programs  to  handle  space  charge 
effects,  they  are  totally  separate  from  those  programs  used  in  lens  or  deflector  design  while  in  reality  you  may  want 


to  include  space  charge  effects  is  the- leas  design  (Ufe.  There  is  no  oae  program  which  completely  models  an 
electron  hem  system  from  source  all  the  way  through  optical  etemema  and  finally  sample  interaction. 

An  area  that  current  computer  modeling  tails  ahott  in  it  optimization  tools.  Ideally,  one  should  give  an 
approximate  leas  design  to  the  rrenputrr  and  the  computer  would  then  carry  out  a  search  through  parameter  space  to 
try  and  find  a  more  optimal  aolubou.  At  present,  such  software  is  limited.  Munro  has  put  out  the  programs  LITHO 
[68]  and  ST1G  [20],  which  performs  optimisation  of  stigmatnr  and  focusing  systems,  but  the  optimization  is  limited 
in  timer  parameters  Him  electrode  pmrotisli  tnd  roil  mnmfi  wtiirh  rlo  not  rrquirr  new  firlrl  calculations  Other 
reaearchew  have  investigated  geometrical  variations  in  the  leases,  but  they  use  many  approximations  to  evaluate  lens 
performance.  Typical  is  van  der  Steen  tt  al  [69],  Szilagyi  and  Szep  [70],  and  Lenz  [71]  who  solve  for  the  axial 
potential  of  a  cylindricaQy  symmetric  leas  using  various  approximations  and  then  use  the  paraxial  ray  equation  for 
calculating  trajectories.  Kato  and  Tsuno  [72]  have  attempted  an  optimization  which  modifies  the  geometrical  shape 
of  the  lens,  but  they  limited  their  optimization  to  minimizing  the  spherical  aberration  coefficient.  Vertes  et  al  [73] 
discuss  lens  optimization  for  analytical  instruments  (e.g.  mass  spectrometry);  other  work  may  be  found  in  Olson  and 
Kusse  [74],  and  Rayces  and  Lebich  [75].  The  most  difficult  obstacle  to  overcome  in  finding  a  global  optimization  is 
the  prodigious  amount  of  computing  power  required.  One  algorithm  which  has  become  popular  for  solving  multi¬ 
dimensional  optimization  problems  is  the  simulated  annealing  algorithm;  Forbes  and  Jones  [76]  discuss  this  method 
as  applied  to  lens  optimization. 

Lastly,  although  automated  mesh  generation  programs  exist  (c.g.  the  program  SOGUN),  they  are  typically 
not  fully  autonomous.  Instead,  the  user  (who  will  have  insight  into  the  problem  to  be  solved)  must  still  instruct  the 
program  where  to  put  the  heaviest  concentration  of  grid  lines.  This  may  still  be  a  tedious  process.  Furthermore, 
there  is  probably  always  a  need  for  more  sophisticated  numerical  methods  to  get  ever  more  accurate  field  calculations. 
One  technique  that  may  relieve  the  user  from  guiding  the  mesh  program  while  also  being  an  accurate  potential  solver 
is  the  use  of  adaptive  grids.  To  illustrate,  suppose  you  allow  the  computer  to  autonomously  set  up  a  coarse, 
uniformly  spaced  grid  and  run  an  FEM  code  once  to  obtain  a  preliminary  answer  for  the  potentials.  The  computer 
could  look  to  find  where  the  equipotential  lines  lie  in  this  first  solution  and  then  completely  readjust  the  grid  so  that 
now  the  grid  lines  conform  to  the  equipotential  lines.  This  would  automatically  solve  the  problem  of  ensuring  that 
enough  grid  elements  are  present  where  the  field  is  rapidly  changing.  A  few  iterations  should  converge  to  an 
extremely  accurate  answer.  Furthermore,  there  may  be  fundamental  numerical  accuracy  benefits  or  new  algorithms 
that  can  be  exploited  if  the  grid  lines  are  on  approximate  equipotentials. 
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